Genome-wide association mapping and genomic prediction of agronomical traits and breeding values in Iranian wheat under rain-fed and well-watered conditions

Background The markers detected by genome-wide association study (GWAS) make it possible to dissect genetic structure and diversity at many loci. This can enable a wheat breeder to reveal and used genomic loci controlling drought tolerance. This study was focused on determining the population structure of Iranian 208 wheat landraces and 90 cultivars via genotyping-by-sequencing (GBS) and also on detecting marker-trait associations (MTAs) by GWAS and genomic prediction (GS) of wheat agronomic traits for drought-tolerance breeding. GWASs were conducted using both the original phenotypes (pGWAS) and estimated breeding values (eGWAS). The bayesian ridge regression (BRR), genomic best linear unbiased prediction (gBLUP), and ridge regression-best linear unbiased prediction (rrBLUP) approaches were used to estimate breeding values and estimate prediction accuracies in genomic selection. Results Population structure analysis using 2,174,975 SNPs revealed four genetically distinct sub-populations from wheat accessions. D-Genome harbored the lowest number of significant marker pairs and the highest linkage disequilibrium (LD), reflecting different evolutionary histories of wheat genomes. From pGWAS, BRR, gBLUP, and rrBLUP, 284, 363, 359 and 295 significant MTAs were found under normal and 195, 365, 362 and 302 under stress conditions, respectively. The gBLUP with the most similarity (80.98 and 71.28% in well-watered and rain-fed environments, correspondingly) with the pGWAS method in the terms of discovered significant SNPs, suggesting the potential of gBLUP in uncovering SNPs. Results from gene ontology revealed that 29 and 30 SNPs in the imputed dataset were located in protein-coding regions for well-watered and rain-fed conditions, respectively. gBLUP model revealed genetic effects better than other models, suggesting a suitable tool for genome selection in wheat. Conclusion We illustrate that Iranian landraces of bread wheat contain novel alleles that are adaptive to drought stress environments. gBLUP model can be helpful for fine mapping and cloning of the relevant QTLs and genes, and for carrying out trait introgression and marker-assisted selection in both normal and drought environments in wheat collections. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-022-08968-w.

2]. Genetic research on this crop has led to its improved productivity. For example, the last decade (2011-2020) witnessed ~ 1% yield increase per annum [3]. However, further improvement is imperative to feed the global population, which will reach over 9 B by 2050 [4]. As the most important detrimental factor, wheat production is restricted by water-limited conditions in most parts of the world. Improvement of crop tolerance to drought stress is one of the essential efforts that can guarantee sustainable yield in wheat fields [2,4]. Right now, research attempts are focusing on exploring the genetic foundation of drought tolerance traits by using association analysis of agronomic characteristics and genomic regions [5].
The breeding of high-yielding and drought-tolerant wheat varieties continues to be a challenging task, because of large "environment×genotype" interactions and low heritability related to yield as a complicated agronomic property [6]. To overcome this problem, highthroughput methods in phenomics, including digital imaging, and in genomics, including association mapping, have been used to uncover the genetic mechanisms underlying yield and its relative characteristics under drought. The findings obtained from these methodologies had been practical for further enhancement in wheat yield not only in water-restricted environmental conditions but also in drought-stressed environments [3].
The advent of next-generation sequencing technologies has provided an opportunity to evaluate genetic variation and discover new markers through implementing the genotyping-by-sequencing (GBS) approach [7]. From this approach, molecular markers such as single nucleotide polymorphism (SNP) have been successfully adopted to discover the complicated agronomical properties of wheat and also have been well-known as key elements in the genome-wide association study (GWAS) approach [8]. The purpose of this approach is to detect genomic regions that can either be QTL, gene, or marker related to important traits for gene introgression, gene discovery, or marker-assisted breeding [2]. The markers detected by GWAS make it possible to dissect genetic structure and diversity at many loci. This can enable a wheat breeder to reveal and used genomic loci controlling drought tolerance [5].
In addition to trait mean-based GWAS (pGWAS), there is a chance to estimate breeding values by some methods such as BRR (bayesian ridge regression), gBLUP (genomic best linear unbiased prediction), and rrBLUP (ridge regression-best linear unbiased prediction) and use them in association mapping (i.e., eGWAS). There is a lack of certainty on the best algorithm when utilizing a multiple-regression model in genomic selection and GWAS since the structure of the population and the architecture of the trait have a remarkable effect on identifying marker impacts [9]. As a result, it is imperative to compare the findings from the various algorithms when dissecting the genetic basis of a complicated trait in a crop population for the first time. This process ensures the efficient detection of QTLs responsible for controlling a quantitative trait, and better control of the error of type I, which is often higher in association mapping studies [10].
To date, about 800 marker-trait associations (MTA) and quantitative trait loci (QTL) have been discovered for wheat drought tolerance traits, including yield, root, physiological, and agronomic ones by using association mapping (~ 100 MTAs) and bi-parental mapping (~ 700 QTLs). Only 70 loci, however, are known as the major genomic regions explaining more than 20% of phenotype diversity [11]. In the past, association mapping research in drought-stressed wheat has utilized a small number of molecular markers [12][13][14][15][16], which seems inadequate for efficiently exploring diversity in diverse wheat collections.
Genomic prediction (GP) is a powerful tool to boost the efficiency and speed of breeding schedules by reducing time cycles and increasing selection accuracy. This approach provides an opportunity by which a candidate gene can be chosen via genotyping before phenotype determination [17]. Genomic prediction utilizes all genetic markers within a model to train a prediction model, which is consisted of all genetic impacts. The model is applied to a validation set for estimating its accuracy [18]. Several studies have demonstrated high or moderate GP accuracy for quantitative characteristics in barley (Hordeum vulgare L.) [19], maize (Zea mays L.) [20], rice (Oryza sativa L.) [21], oat (Avena sativa L.) [22] and wheat (Triticum aestivum L.) [17].
This study was aimed at detecting drought tolerance candidate QTLs, genes, or markers linked with agronomical traits by using GWAS in 208 wheat landraces and 90 cultivars grown under normal and drought conditions. In eGWAS, the goal is to identify SNPs related to the correction value of the traits, which are passed on to the next generation. The next purpose of this work was to select the best model for estimating prediction accuracies in genomic selection. To the best of our knowledge, our report is the first study on pGWAS and eGWAS of agronomical characteristics in Iranian wheat landraces under rain-fed and wellwatered conditions. The findings from this research will be an interesting source for marker-assisted breeding, genomic selection, introgression of favorable genes into high-yielding cultivars, and improvement of yield-associated characteristics under drought.

Phenotypic data summary
In this study, 298 landraces and cultivars of bread wheat were grown under rain-fed and well-watered conditions and analyzed for various agronomic traits. According to the analysis of variance, genotypic, environmental, and genotype×environmental effects on agronomical traits were significant under rain-fed and well-watered environments. Variances associated with genotypic effects were higher than those associated with environment and genotype×environment effects across all traits, indicating genotypic effects had a greater impact. There is a high heritability in plant height traits, but a low heritability in grain yield traits. However, the agronomical traits of wheat grain showed acceptable heritability (Table S1). The box plots related to eight agronomical traits of wheat landraces and cultivars under favorable conditions (well-watered) and drought stress (rainfed) are shown in Fig. 1. The mean of all traits under stress decreased when compared to a normal situation in both cultivars and native populations implying the presence of considerable diversity in agronomical traits of wheat accessions, and this variation is greater in native populations. The mean of all traits, except plant height, in both conditions, was higher in cultivars than in landraces.

Clustering analysis
Under normal conditions, the heatmap was plotted based on the mean of agronomic traits and breeding values by using three methods: BRR, gBLUP, and rrBLUP. From the results, wheat accessions were clustered into four groups. In clustering based on the mean of traits, Group No.1 included 82 high-yielding genotypes that were 41 cultivars and 41 landraces, Group No.2 consisted of 89 genotypes with average to high yield (24 cultivars and 65 landraces), Group No.3 contained 44 genotypes with average to low yield (21 cultivars and 23 landraces), and Group No.4 composed of 83 low yielding  genotypes that were mainly native populations (4 cultivars and 79 landraces) (Fig. 2a). From the BRR method, wheat genotypes were divided into four groups; the first, second, third, and fourth groups consisted of 61, 42, 104, and 91 genotypes, respectively (Fig. 2b). From the gBLUP, the first group included 85 genotypes with a high breeding value of grain yield (72 cultivars and 13 landraces), the second group consisted of 102 genotypes with medium to high breeding value for yield and yield components (16 cultivars and 86 landraces), the third group contained 97 genotypes with medium to low breeding value for yield and components (2 cultivars and 97 landraces), the fourth group composed of genotypes (17 landraces) with low breeding values for yield and yield components (Fig. 2c). From the BRR method, wheat genotypes were divided into four groups; the first, second, third, and fourth groups consisted of 69 (67 cultivars and 2 landraces), 59 (9 cultivars and 50 landraces), 88 (12 cultivars and 76 landraces), and 82 genotypes (2 cultivars and 88 landraces), respectively (Fig. 2d). The results of gBLUP were most similar to the trait mean method in terms of genotype clustering.
Drought-stressed genotypes were also classified into four groups based on the trait mean and the breeding value methods. In clustering based on the mean of traits, the cluster 1 included 31 genotypes with high yield, which were mainly cultivars (18 cultivars and 13 landraces), the cluster 2 consisted of 123 genotypes with average to high yield (24 cultivars and 99 landraces), the cluster 3 contained 43 genotypes with average to low yield (19 cultivars and 24 landraces), and cluster 4 composed of 101 genotypes with low average yield, which were mainly native populations (29 cultivars and 72 landraces) (Fig. 3a). From the BRR, the first group included 61 cultivars with a high breeding value of grain yield, the second group consisted of 67 genotypes (18 cultivars and 49 landraces) with medium to high breeding value for yield and yield components, the third group contained 53 genotypes with medium to low breeding value for yield and components (8 cultivars and 45 landraces), the fourth group composed of 117 genotypes (3 cultivars and 114 landraces) with low breeding values for yield and yield components (Fig. 3b). From the gBLUP method, wheat genotypes were divided into four groups; the first, second, third, and fourth groups consisted of 65, 83, 48, and 102 genotypes, respectively (Fig. 3c). Clustering based on breeding values by using BRR, gBLUP, and rrBLUP had 42, 48, and 39% similarity in terms of genotype clustering in different clusters, respectively. This indicates that the gBLUP categorized wheat accessions more accurately than the other two BRR and rrBLUP methods (Fig. 3b, c and d).

Linkage Disequilibrium (LD)
LD assessment indicated that this indicator varies between chromosomes and across each chromosome and it usually decreases with increasing distances between SNP locations. A total of 1,858,425 marker pairs with r 2 = 0.211 were identified in cultivars, of which 700,991 (37.72%) harbored significant linkages at P < 0.001. The strongest LD was recorded between marker pairs on chr 4 A (r 2 = 0.318). Genomes D and B possessed the lowest (63,924) and highest (370,359) number of significant marker pairs, respectively. A similar assessment on wheat landraces found 1,867,575 marker pairs with r 2 = 0.182, of which 847,725 (45.39%) harbored significant linkages at P < 0.01. Similar to cultivars, marker pairs on chr 4 A showed the strongest LD (r 2 = 0.369). Genomes D and B possessed the lowest and highest number of marker pairs (92,702 and 427,017), respectively. In the D genome, the LD decay was slower than the LD decay in A and B genomes, indicating that the size of the linkage blocks is larger in the D genome. In addition, in cultivars, compared to the native populations in genome D, the LD decay was slower, which probably indicates the selection of more genome-related traits in breeding work. Based on the observations, the most significant marker pairs in wheat landraces were found at distance < 10 cM (Table 1).

Population structure
To estimate the subpopulations, the ΔK value was plotted against the number of clusters (K). The largest ΔK value was found at K = 3, reflecting three population substructures, Sub.1, Sub.2, and Sub.3 (Fig. 4a). Sub.1 included 113 genotypes with 6 cultivars and 107 landraces; Sub2 contained 111 genotypes with 97 landraces and 14 cultivars; Sub.3 consisted of 74 genotypes with 70 cultivars and 4 landraces (Fig. 4b). From PCA analysis, the estimated PCs showed that PCs 1 and 2 explained 10.29 and 6.28% of the genotypic variation, respectively (Fig. 4c). Cluster analysis using the kinship matrix also supported the STRU CTU RE results (Fig. 4d).

Genome-wide association studies for agronomic traits and estimated breeding values
Under optimal irrigation and using imputed markers and -log10 P > 3, 283 significant SNPs were discovered for agronomic characteristics by MLM. Of these, 106, 137, and 40 markers were for genomes A, B, and D, respectively. Therefore, genome B had the highest number of significant SNPs. The number of significant markers for PH, GY, GN, TKW, SW, SA, SH, and SF were 39,57,19,48,11,31,43, and 35, respectively (Fig. S2a). The number of significant SNPs based on BRR, gBLUP, and rrBLUP were 362, 358, and 294, respectively. (Fig. S2b, c and d) The gBLUP method with the most similarity (81.27%) in  the terms of significant markers had the best justification when compared to other methods (Table 2)  In stress, less significant markers were identified than the normal situation, 194 significant SNPs were identified by the MLM method; Of these, 48, 129, and 17 markers belonged to genomes A, B, and D, respectively. Genome B had the highest percentage of significant SNPs in a stressful environment. The number of significant markers for PH, GY, GN, TKW, SW, SA, SH, and SF were 9, 30, 16, 21, 15, 31, 31, and 41, respectively (Fig. S4a). The number of significant SNPs obtained by BRR, gBLUP, and rrBLUP methods was 364, 361, and 301, respectively (Fig. S4b, c and d). The gBLUP with the most similarity (71.64%) in the terms of significant markers had the best justification when compared to other methods (Table 2). By BRR, gBLUP, and rrBLUP, a total of 134, 121, and 97 significant SNPs for genome A, 187, 198, and 167 SNPs for genome B, as well as 43, 42, and 37 SNPs for genome D were identified, respectively (Fig. S4b, c and d). The Manhattan circular plot shows significant SNPs at P value < 0.001 (black) and < 0.00001 (red) (Fig. 6). The Manhattan rectangular and Q-Q plot are shown in Fig. S5.

Gene ontology
The markers with the highest significance (P < 0.0001) and pleiotropic impact were studied in more detail. In the normal environment, 29 markers containing overlapping genes were identified that are involved in important biological and molecular processes. 12 markers were identified based on the pGWAS method and 17 markers were identified based on the eGWAS method. The number of GO based on BRR, gBLUP, and rrBLUP were 18, 15, and 16, respectively. The gBLUP and BRR method was most similar to (66.67%) the pGWAS method. The most significant markers were located on chr 6B, 5B, and 5 A. Of these, 8 SNPs were detected by both pGWAS and eGWAS methods. Some of the uncovered MTAs were responsible for the following molecular and biological processes: lipid biosynthetic process, protein-binding, carbohydrate-binding, lipid transport, RNA-binding, protein ubiquitination, protein deubiquitination, protein catabolic regulation, nucleoside metabolic process, UMP salvage, CTP salvage, and ubiquitin-dependent protein catabolic process (Table 3).

Genomic prediction
The gBLUP, rrBLUP and BRR approaches using imputed SNPs led to the identification of the highest prediction accuracies for 5, 3, and 1 phenotypes in rain-fed, and 5, 3, and zero phenotypes in well-irrigated environments, respectively (Fig. 7). Under rain-fed, the highest prediction accuracy was determined via the gBLUP model for GY  (Fig. 7).

Discussion
Shedding light on the genetic mechanisms controlling quantitative traits such as grain yield in wheat represents an opportunity for the improvement of drought tolerance. To achieve this goal, this experiment aimed at exploring the structure of the population and at uncovering MTAs in Iranian wheat accessions. Significant, positive correlations among the wheat characteristics confirmed the value of the data in the current GWAS analysis. This is evidenced by Laido et al. [26] who highlighted the relationship between morphological characteristics having a high correlation to detect relevant QTLs.       High correlation occurring between agronomic traits can be justified by indirect or direct contributions of one trait to another [27]. Taking a look at the wheat genome, genomic regions responsible for such agronomic characteristics can be equivalent. This is supported by the presence of multi-trait correlations where one gene has a pleiotropic impact on highly-associated characteristics [2]. For example, Mwadzingeni et al. [8] showed that one locus controls several wheat properties such as grains per spike, spike length, and plant height, which are highly linked often [28]. Such observations support the requirement to confirm if such locus is not also linked to another trait, because it shares similar sequences with the regions responsible for the latter trait. Some loci, however, affect only one crop property [8].
Breeding value-clustering by using BRR, gBLUP, and rrBLUP had 77, 68, and 83% similarity with the trait mean method in the terms of wheat accessions grouping, respectively. This indicates that rrBLUP can categorize wheat accessions more accurately than the other methods. Moreover, rrBLUP with the most similarity with the trait mean method in the terms of discovered significant markers, suggesting its potential in uncovering SNPs. As a result, rrBLUP model can detect genetic impacts in wheat populations better than other models.
Overall, obtaining the best outcomes from the breeding value-based methods depend on the genetic architecture of trait, genetic variation, etc. [18].

Linkage disequilibrium of markers
Of the results, the SNPs covered the wheat genome well. The SNPs were higher in genome B. The higher frequency of SNPs in genome B results from the evolutionary events [29]. Genomes D had the highest LD followed by genome A, followed by genome B. At the chromosome level, the strongest LD was recorded between marker pairs on chr 4 A. The fact that cultivars exhibited higher LD in contrast to landraces, particularly in the genome D, is presumably a consequence of selection throughout the time of breeding efforts [30]. The presence of closely linked marker pairs with non-significant LDs and marker pairs in LD over a long distance in this research has been shown previously in wheat and other crops [8,31]. This reflects that LD is not static because LD can be affected by various elements including genetic admixture [8].

Population structure of Iranian wheat accessions
The population under consideration was divided into four distinct sub-populations. This is expected because the wheat accessions have diverse pedigrees. Of course, the presence of common parents or origins in the pedigree of accessions often leads to some relationships among them [2]. The findings derived from the population substructure analysis are beneficial in following superior parents that can be used in the improvement of wheat tolerance to drought stress conditions [3]. Therefore, latter researchers can utilize this genetic pool to employ the genetically disparate accessions, which in turn exhibit wheat farmer-preferred properties.

SNPs and MTAs for wheat agronomic traits
From a brief look at the number of SNPs, lower significant SNPs were recorded under drought than normal conditions, reflecting GWAS analysis for exploring drought tolerance is affected greatly by environment*genotype interactions [8].
This experiment led to discovering of a total of 29 and 30 highly significant MTAs in normal and drought environmental conditions, correspondingly. Albeit only those associations at P < 0.0001 were regarded as significant, the rest of these MTAs may be helpful for enhancing wheat tolerance to drought stress. These associations can be located in genomic regions affecting the agronomic characteristics. The MTAs for yield appeared significant at a higher P value, because this trait is highly complicated in genetic nature with low heritability [32].
To date, many attempts have been focused on locating QTLs and genes affecting wheat traits in drought environments for facilitating marker-assisted breeding [2, 3]. The MTAs detected in this study are added to the previous pool of candidate genes and markers. However, it is a challenging task to align our results with earlier works because of the use of disparate reference genomes than the IWGSC Ref.Seq, the lack of accurate genomic locations, or the utilization of various markers (GBS-derived SNP vs. SSR and DART) [2,3,5,9]. Of course, detection of MTAs on the same chromosome as previous projects increases the assurance of these MTAs.
Four MTAs for grain yield were recorded on chr 3B, 4 A, 5 A, and 3D in this study. Earlier research efforts have discovered MTAs/QTLs for grain yield on wheat chr 7B [31,33,34], 7 A [31,[34][35][36], 5B [15,31,34], 3D [34], 3 A [31,34,37,38], 2B [34,[37][38][39][40], and 1B [34,38,39]. Thus, MTAs on chr 3B, 4 A, and 5 A have not been reported and they are new for wheat yield. Six MTAs for TKW were found on chr 5 A, 1B, 3B, 6B, 1D, and 2D. Earlier reports have detected MTAs/QTLs for TKW on chr 7D [35], 7B [31], 5B [41], 3B [35], 3 A [40,41], 2D [39], 2B [31,35,39,42], 2 A [35], 1 A [31,[39][40][41] and 1B [43]. For plant height, two MTAs were revealed on each of chr 5B, 6B, and 2D. All 21 chromosomes carry genes that control plant height in wheat [42,44,45]. Up to now, 24 reduced height (Rht) genes (Rht1-Rht24) are catalogued in wheat [46,47], where Rht8 on chromosome arm 2DS has been extensively explored [48,49]. We could locate only two QTLs to chromosome 2DL, whereas the ones reported by Borner et al. [50], on chromosome 2DS could not be detected. Other MTAs detected in our research effort were responsible for grains per spike, spike weight, spike fertility spike area, and spike harvest index. Some of the MTAs detected in this study were involved in the following important biological and molecular processes: metal ion binding, monooxygenase, acyltransferase, oxidoreductase , acyltransferase, methyltransferase, peptidase, and dependent microtubule motor activity. The gBLUP with the most similarity (80.98 and 71.28% in wellwatered and rain-fed environments) with the trait mean method in the terms of discovered significant SNPs, suggesting the potential of gBLUP in uncovering SNPs. The results show that the gBLUP method performs better than the rrBLUP and BRR methods in terms of predicting the accuracy of genomic breeding values. In gBLUP, genomic relationships are used to estimate an individual's genetic merit. Genomic relationships are estimated based on DNA marker information for this purpose. To make better predictions of merit, the matrix defines the covariance between individuals on the basis of observed similarity rather than expected similarity based on pedigree. Several studies have described the gBLUP method for estimating genomic breeding values [51][52][53][54]. Research shows that gBLUP and rrBLUP are similar models. One of the advantages of gBLUP over rrBLUP is the reduction of the dimensions of the mixed equations to the number of people in the reference population, the calculation of accuracy and error predicting corrective values as commonly used in pedigree methods and combining The information of genotyped and non-genotyped individuals was mentioned simultaneously in the mixed equations [18].
Based on the GO results, the BRR and gBLUP methods were able to better identify the relationship between the studied traits, respectively, and were most similar to the pGWAS results. Generally speaking, genes/markers affecting a trait under drought also are responsible for that trait under normal conditions [8]. Ideally, the impacts of such genes/markers may not be influenced by any moderate changes in environmental conditions, thus they can be helpful in gene introgression or markerassisted selection when adaptation improvement [55]. Some genes/markers, on the other hand, may affect specific traits differentially under various conditions [55].
Our findings suggested that genomic prediction is a helpful tool for predictive characterization of wheat genotypes, permitting phenotyping to be limited to a fraction of the germplasm rather than the whole collection [56][57][58]. Similarly, Kehel et al. [59] stated that genomic selection can be used within wheat accessions to predict key traits with an accuracy of more than 0.7, more especially for the traits with high to moderate heritability. Accounting for stratified populations is usually carried out by the first five principal components as covariates in a prediction model [57,60,61]. As expected, a significant population structure was identified in the Iranian wheat landraces, with the first five eigenvalues accounting for 30.5% of genetic diversity. The population structure indicated a negative effect on performance in GWAS and GP models, which was also exhibited in other researches [61,62]. Of our observations, the highest prediction accuracy was achieved via the gBLUP model. Shabannejad et al. [18] evaluated classic approaches for exploiting GP accuracy by BRR, gBLUP, rrBLUP models in normal and drought environments in wheat cultivars and landraces. They identified the highest GP accuracies via the gBLUP and BRR method. The authors observed that obtaining the highest GP accuracy depends on the genetic variation, genetic architecture of trait, level of LD, and the genomic selection approach. As a result, the gBLUP model can detect genetic impacts in wheat populations better than other genomic prediction models.

Conclusion
MTAs are the key elements to detecting genomic regions related to wheat agronomic traits under drought stress. The current experiment found 29 and 30 highly significant MTAs under normal and drought conditions. The markers detected would be useful genomic sources for cloning and fine mapping of underlying genes, and for conducting gene introgression and marker-based selection in wheat under normal and drought conditions. A further research attempt is needed for validating the markers detected in the current project using a larger wheat population.

Plant material and experimental conditions
A field research effort was performed in two growing seasons (2018-19 and 2019-20) under rain-fed (drought) and well-watered (normal) conditions at the research farm, University of Tehran, Iran. In this study, 90 cultivars and 208 landraces (Table S2) of wheat were investigated in an alpha-lattice experiment with two replications. The wheat accessions were cultivated in the plots including four rows (1*1 m 2 ) at 0.5 m intervals. In the well-watered crops, the threshold of irrigation was regarded based on 40 mm evaporation from a standard pan. The reference crop evapotranspiration [ET 0 = E pan × K pan ; where K pan is a pan coefficient (0.8) for each month and E pan is the evaporation depth from the pan surface (40 mm)] and crop coefficient [K C ] were estimated to measure evapotranspiration (ET C = K C × ET 0 ) [63]. The time of irrigation was determined from the ratio of the assigned water for 1400 m 2 (the cultivation area of total genotypes in two replications) to water discharge (10.8 m 3 /h). The volume of water required for each hectare (m 3 /ha) was calculated via the depth of ET 0 (mm) multiplied by ten. The rain-fed crops were exposed to rainfall, which was the only accessible water source. The monthly rainfall pattern for the growing seasons is represented in Table  S3. At the maturity stage, 20 plants were harvested from the middle rows of plots to measure traits, including spike fertility (ratio of grain number to spike weight), thousand-kernel weight (g), grain yield (g per plant), grain number per spike, spike weight (g), spike harvest index (ratio of spike grain weight to spike weight, %), spike area (cm 2 ), and plant height (cm).

GBS analysis
To sequence wheat accessions, this experiment followed the procedure as explained by Alipour et al. [29] to establish the GBS libraries. After trimming reads to 64 bp and categorizing them, single nucleotide polymorphisms were discovered by internal alignment. SNPs were called through the UNEAK GBS pipeline, where SNPs with lowallele frequency < 1% and low-quality scores < 15 were discarded to reduce false positives. The SNP imputation process was implemented by available allele frequencies in BEAGLE V.3.3.2 [64]. The LD was calculated by the TAS-SEL V.5 [65]. The W7984 reference genome was adopted in the recent study because of fulfilling the highest accuracy of imputation among the wheat references [30].

Structure of wheat population
Population structure in the Iranian wheat accessions was revealed by STRU CTU RE V.2.3.4. In this software, the parameters were set at 30,000 burn-in periods, with 30,000 MCMC iterations after burn-in [66]. To permit the picking up of repetition with the highest value of Ln likelihood, 10 replications were run for K values of 1 to 10. By using TASSEL software, genotypic data of wheat accessions were imputed [67]. Moreover, principal component analysis (PCA) was conducted to verify the STRU CTU RE outcome. To determine the accession relationships, a neighbor-joining analysis was carried out by TASSEL V.5. Linkage disequilibrium (LD) was determined through R 2 value, squared allele frequency correlation, from which the significant allele pairs were estimated by 1,000 permutations.

Trait mean-based GWAS (pGWAS)
The mixed linear model (MLM) was followed to estimate the marker impacts on the wheat population. The general linear model was conducted by population structure matrix (Q) integrated as a covariate for correcting the effect of subpopulations. The mixed linear model was performed by both the family structure matrix (Kinship, K) and Q for controlling both errors of type I and II. The association mapping was implemented using MLM functions of TASSEL V.5. To correct for multiple test, a false discovery rate was utilized to declare significant MTAs [66,68]. For a better answer in the recent study, only the outcomes of the MLM procedure were given. There are several methods to determine the threshold in GWAS and all of them have some advantage and disadvantage. But, the most important thing is confirming the results using further analysis. Here the threshold -logP > 3 was considered to find higher number of significant SNPs and identify the important ones using GO and pathway analysis. While from the threshold of -logP > 5 was considered to identify very significant and important SNPs. To explore associations between genotype and phenotype, a Manhattan plot was obtained using the CMplot package [69].

Breeding value-based GWAS (eGWAS)
Three methods rrBLUP [70], BRR [71], and gBLUP [72] using the Intelligent Prediction and Association Tool (iPat) software were used to obtain the breeding values. A mixed linear model (MLM) was used to estimate the effects of markers using breeding values on wheat populations [9].

Annotation of putative candidate MTAs
The ensemble-gramene database was employed to extract the molecular and biological functions of SNPs in the gene ontology by using the IWGSC RefSeq V.2.0, which has been provided for the Chinese Spring [http:// www. grame ne. org/]. Furthermore, the significant SNPs were analyzed via KOBAS version 2.0 for gene ontology enrichment analysis in KEGG [https:// www. genome. jp/ kegg/].

Genomic prediction strategies
GP was calculated by various approaches: BRR [71,73], gBLUP [72,73], and rrBLUP [70,73]. All of the analyses were performed by iPat [74]. For the population, 20% of genotypes were assigned randomly to a validation set and all of the residuals were utilized as a training set. This process was reiterated 100 times for all of the prediction approaches. The GP accuracy was calculated as Pearson's correlation (r) between BLUPs and GEBVs over the validation and training sets [75].

Statistical analysis
The descriptive statistics and correlation analysis were implemented by R V.4.1 using the dplyr, ggpubr, psych, and ggplot2 packages. Heatmap analysis was carried out using heatmap.2 function in gplots R package to classify wheat accessions.